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ABSTRACT 

We demonstrate the high accuracy of the density spHtting method to compute the gravita- 
tional potential and field in the plane of razor-thin, axially symmetric discs, as preliminarily 
outlined in Pierens & Hure (2004). Because residual kernels in Poisson integrals are not C°°- 
class functions, we use a dynamical space mapping in order to increase the efficiency of advanced 
quadrature schemes. In terms of accuracy, results are better by orders of magnitude than for the 
classical FFT-methods. 

Subject headings: gravitation — methods: numerical 

1. Introduction and motivations 

Discs are ubiquitous objects in the Universe and span different velocity/length scales: galactic (stellar) 
discs. Active Galactic Nuclei discs, circumstellar discs, binary and circum-binary discs, sub-nebulae. For 
many of them, self-gravity plays role in their structure and dynamics and so, the gravitational potential 
and associated accelerations are required at a certain level of disc modeling. Solving the Poisson equation 
in extended, continuous media like gas discs is however not trivial practically, and it has occupied astro- 
physicists for many decades. For time-dependent simulations, fast but low accuracy algorithms are generally 
preferred. For steady state analysis, accuracy is more critical than computing time, as it allows for instance 
to characterize the precise connection between various branches of solutions (e.g. Hachisu 1986, Ansorg, 
Kleinwachter & Meinel 2003). 

Many numerical methods have been proposed, but a very few uses the integral formalism. In a previous 
paper (Pierens & Hure 2004; hereafter Paper I), we have outlined a method to avoid the singularity in the 
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Poisson kernel so that the field in the plane of razor-thin, axially symmetric discs can easily be accurately 
computed from elliptic integrals by a single radial quadrature. The motivation of the present paper is 
twofold. First, Paper I just touches the problem at the theoretical level without giving numerical examples 
and discussing the implementation and possible performances of the method. Also, the potential was not 
considered. This is done here. Second, we have recently realized that the classical FFT-method (e.g. Binney 
& Tremaine 1987) which is among the most widespread methods, have a much lower precision and lower 
order of convergence, in comparison. We show actually here that the density splitting method can exhibit a 
high convergence order, depending on the quadrature rule, and the precision can easily reach the machine 
precision with a few tens sources points, at least for smooth surface density profiles. 

We outline the splitting method for the gravitational potential and radial field in Section 2. We then 
stress the non-derivability of Poisson kernels in Section 3, and propose a space mapping. In Section 4, we 
illustrate the possible performances of the method on a test-case (namely, a finite size disc with exponentially 
decreasing surface density profile), with three different quadrature rules. We compare the accuracy of the 
method with the classical FFT-method in Section 5. A few concluding remarks are found in the last section. 




2. Outline of the density splitting method 

The gravitational potential \1/ and radial component of the field gj^ = — 9/?^' in the plane of a razor thin 
axially symmetric disc (see Fig. 1) is""^ (e.g. Durand 1964): 




where is the inner edge, Uout is the outer edge, S(a) is the surface density, 

are the Poisson kernels, K and E denote the complete elliptic integral of the first and second kinds respec- 
tively, k = 2\J aR/{a -f i?) < 1 is their modulus, and vj — {a — R)/{a + R). As it is well known, these 
kernels diverges when a R, with the result that, in practice, neither ^' nor gn can properly be determined 
by direct integration for R e [ain,aout]- As outlined^ in Pierens & Hure (2004) (hereafter Paper I), such a 



^Matrix notation is employed only for compactness. 
^Only the field was considered in Paper I. 
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singular behavior can be avoided if the surface density profile is split into two components according to 

S(a) = So + (5S(a,i?), (3) 

where Sq = is the local value, and is the remainder (a function which depends on a and R). 

Consequently, the potential and the radial field in the disc are given by 



\]> \ /\J,homo.\ 



,res. 



where vj/homo- a,nd 5^^°™°' are analytical functions (proportional to So; see Appendix A in Paper f for the 
field, and Appendix A in this paper for the potential). Terms and g'^^' correspond to the departure 
from the homogeneous disc. These are simply given by 



\9TJ J a 



""^ <5E ( ') da. (5) 



The point is that both 6T, x Kg and i5S x are finite when a = R, although ^ oo and Kg ^ oo 
(see Appendix B in Paper I for a proof), with 

((52 X K*)|„^^ = 0, 

(6) 

(^Sx^,)U^==2G(f)^^^. 

Note that, in this procedure, the radial derivative of the surface density ^ is needed at each field point when 
computing g^^^- from Eq.(5). This is quite uncomfortable if the surface density is not defined analytically, 
but on a grid as in many simulations. 



3. Residual integrands are not C°°-class functions 

Acciiracy of residual terms is set by the scheme performing the numerical quadrature, provided the 
integrand is a well-behaved function, fn the present problem, a difficulty arises because the residual inte- 
grands (5S X Kg and 5Y, x k-h are continuous, but not C°°-class functions (i.e. differentiable for all degrees 
of differentiation). Even, these are non-derivable (i.e. not C°). This is easily understood from the first 
derivative 

d , _^ dK 

((5S X k) = — X «;-h(5E X — , (7) 



where k denotes either or «;„ 



and 



da da da 



dK dK dk 
da dk da 



- -— (9) 
a (a + i?)2 2a ' ^ ' 




We see that the first term in the right-hand-side of Eq.(7) brings a diverging contribution for a = i? 
since dT,/da can not be zero on the whole integration range. We thus have ^ {6T, x k) ^ oo at the field 
point. This point is illustrated in Fig. 2 which displays the function 6T, x k for both the potential and field 
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Fig. 2. — Top: the function i5E x k$ versus a {dashed line) and its first derivative {plain line) in a disc with 
inner edge ain = 0, outer edge Oout = 1, and exponentially decreasing surface density profile. The field point 
is located at _R = 0.5. Bottom:, same but for the radial field. Areas under ST, x and ST, x Kg {filled zones) 
are the residual terms \E''''"^- and g'^^^- respectively. 



as well as their first derivative in a disc with exponentially decreasing surface density profile. Although not 
visible at the scale of the graphs, there is a small "knee" just at the field point where the first derivatives are 
infinite. Further, because successive derivatives of elliptic integrals inevitably produce the K-function which 
is logarithmically diverging as a ^ i?, we conclude that 

oo for any n > 1. (10) 

a=R 

Values of the integrands SYl x k and its a-derivatives at a = i? are summarized in Tab. 1. It means 
that, if residuals terms are numerically determined following Eq.(5) (i.e. by integration in the a-space), then 
most classical quadrature rules which are based on a certain fitting of the integrand by polynomials cannot 
behave in an optimal manner. This problem especially concerns high-order quadrature rules. We have no 



— {ST, X k) 



Table 1. Values of the integrands x k and its a-derivatives at the field point. 



order of derivative 


JSk^ at a = 




a.t a = R 


0th (function) 





2G{ 


da J a=R 


1st derivative 


oo 




OO 


> 2nd 


oo 




oo 
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idea how to totally remove such a difficulty. It is possible however to reduce it partially by re-sampling the 
integrand in an appropriate way. For instance, if we consider the function w{a) defined by 



^ < if a<R, 

Oout 



w{a) = {O a R = a, (11) 

> if a> R, 



V 



then Eq.(5) now reads 



2aout / ' ("^ ) \w\dw, (12) 



9r^' / J Win V i^g 



where Win = w(ain) and Wout = ''^aout)- The "new integrands" (5E x and ST, x Kg\w\ arc obviously 

regarded as a function of the new variable w. The advantage of this space mapping is twofold. First, it makes 
the first derivative finite; tw-derivatives of the new integrands are summarized in Tab. 2 (see Appendix B for 
a detailed calculus). Figure 3 displays the first, second and third derivatives, in the same conditions as for 
Fig. 2. Second, ST, x Kg\w\ vanishes at the field point (since w = 0). Hence, dT/da is not more needed. 



4. Example of performance 



We briefly show the possible performances of the method through a typical example, namely a disc with 
[oin, Clout] = [0, 1] and an exponentially decreasing surface density profile, as already considered above. Tests 

have concerned a large amount of disc models (various surface density profile, various axis ratios flout /ain)- 
For simplicity, we shall discuss only the potential; results for the field are similar. We consider three different 
schemes to determine at various radii R inside the disc from Eq.(12), namely: 



• the composite trapezoidal rule (hereafter, the CT-rule), with N source points. This is a 2nd-order 
accurate scheme (e.g. Press et al. 1992). 

• a 6th-order, regular-spacing quadrature rule due to Gill & Miller (1972), with A'' source points (here- 
after, the GM-rule). 

• the Gauss-Chebycheff-Lobatto approximation (hereafter, the GCL-rule), with N Chebycheff polyno- 
mials. This collocation method would be the most efflcient technique (with a spectral convergence) in 
the presence of a C°°-class integrand (e.g. Boyd 2001). 



Table 2. Same legend as for Tab. 1, but when in the w-space. 



order of derivative at w = ^SKgluij at w = 



0th (function) 







1st derivative 





T4G(f)„^, 


2nd derivative 





oo 


3rd derivative 


oo 


oo 


> 4th 


oo 


oo 
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Fig. 3. — Same legend as for Fig. 2 but the integrands (5S x x ^ (top) and 6T. x Kg x ^ (bottom). 



We measure the relative precision A\['/^^' on potential values ^ with the e-parameter defined as 



= los 



rcf. 



(13) 



where vp'^^^- is the reference value. Since the exact potential is not known for the case considered, reference 
values \I>'''=f- are obtained by considering a larger number of source points, as commonly done (e.g. Cohl & 
Tohline 1999). We use double precision (DP) computer calculus so that the precision is limited to '--^ 2 x 10^^^ 
(that is e > edp — —15.7). 

Figure 4 gives the results of the splitting method, namely (e) versus N for the three quadrature rules 
listed above, where (e) is an averaged value, obtained by computing ^E* at 2^ equally spaced positions spanning 
the range [ain,aout] (the actual number of field point is unimportant). We see that the accuracy on the 
potential reaches the computer precision for less than a thousand source point with the GM-rule, and for 
only a hundred spectral coefficients for the GCL-rule. Note that the relative precision is better than 0.1% 
with the CT-rule for a dozen source points only, which is remarkable. Also, we observe that the GCL method 
does not exhibit a spectral convergence, for reasons mentioned above but appears very efficient. Indeed, these 
performances are noticeably reduced if quadratures are performed in the a-space. 



5. Comparison with the classical FFT-methods 

It is interesting to compare the accuracy of the density splitting method with that of the classical FFT 
method (e.g Binney & Tremaine 1987) based on N x N polar cells. Results obtained on the same test-model 
are reported on Fig. 4 as open and filled circles (see below). We see that the precision is rather poor even 



-- this work, CT-rule 
— this work, GM-rule 
this work, GCL-rule 
• FFT-method (no smoothing length) 
o FFT-method, smoothing length Ji,=0.01 




Fig. 4. — Averaged e-parameter versus N for the splitting methods and for the FFT-method with smoothing 
length {open circles; A = 0.01) and without smoothing length {filled circles). For the FFT-method, the disc 
inner edge is set to 10~^. See Fig. 5 for the e{R) and A'' = 128. 



for large N, and is definitively lower than that of the splitting method with the CT-rule. 

Let us remind that the FFT-method is probably the most widely method used in simulations of self- 
gravitating discs to compute the gravitational potential (Binney & Tremaine 1987). The only advantage of 
the FFT-method seems to be its great rapidity, N log( A) order in time. Low computing time is a fundamental 
requirement if the Poisson equation is to be coupled with other equations, as it is the case in general. One 
must however realize that the FFT-method has a few major drawbacks. First, it is first-order accurate, 
as can be seen in Fig. 4, that is, one order less than the splitting method with the CT-rule. This means 
that a huge amount of source/grid points is required before reaching great accuracy. If we extrapolate data 
shown in the plot (and ignoring loss of significance and round-off errors that would impose a saturation of 
(e) well above ecp), we find that the FFT-method would need N ~ 10^^ to reach the computer precision 
(that is, 10^^ polar cells). Second, it is a particle-type method (e.g. Hockncy & Eastwood 1988): each cell 
is made homogeneous, and converted into a particle with arbitrary assignment of both location and mass 
density. The precision of the FFT-method is apparent and fortunate (assignment errors almost cancel or 
compensate). Third, the use of the FFT-method in cylindrical coordinates requires a logarithmic spacing 
of grid points. This means that i) the origin can not be included, ii) the outer disc has always a much 
lower resolution than the inner disc (by a factor N), and iii) many interpolations are necessary when other 
equations are solved on a different grid (one has to pass from one grid to the other). It is true that matter in 
astrophysical discs is generally concentrated at the inner edge, but disc self-gravity concerns regions rather 
located near the outer edge. Fig. 5 shows e versus i? for TV = 128 in the example considered before. We 
see that the accuracy is better in the inner disc than in the outer disc with the FFT-method, whereas it is 
uniform with the density splitting method (whatever the quadrature rule). 
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Fig. 5. — e-parameter versus the radius R with the sphtting methods (lines) and with the FFT-method 
{filled circles), for N = 128. For the FFT-method, the disc inner edge is set to 10~^. 

Fourth, the local contribution to the gravitational potential (that is the contribution of the cell containing 
the field point i?) is treated in an artificial way. In this purpose, a smoothing length A is often introduced, 
and its value is chosen arbitrarily. Results obtained for A = 0.01 are shown in Fig. 4 (open circles). Hence, 
potential values are smoothing length-dependent. A slightly better accuracy is obtained with the prescription 
by Binney & Trcmaine (1987) (filled circles in Fig. 4), where the local contribution is estimated analytically, 
without any A. However, as some authors have noticed (e.g. Caunt & Tagger 2001), the FFT-method can 
trigger numerical instabilities in hydrodynamic codes, probably because the order of this method is lower 
(and errors larger) than the order (generally two) of difference schemes used in other fluid equations. 

6. Concluding remarks 

In this paper, we have discussed the possible implementation and performances of the splitting method 
for razor-thin, axially symmetric discs. In particular, we have emphasized the role of a space mapping in 
order to rise the accuracy (or efl'ective order) of advanced quadrature rules. We have noticed the important 
weaknesses of the classical FFT-methods, especially in terms of accuracy, by direct comparisons. Obviously, 
our method is characterized by much longer computing times. Note however that a few grid points are 
necessary to reach a high precision. Besides, the method is easily parallclizable. 

Several extensions and improvements to the present method could be brought. For instance, other 
space mapping than proposed can probably work well. However, it is hard to find a sampling function that 
makes the integrands and derivatives finite (and possibly zero) at the field point without increasing its wings, 
which is numerically uncomfortable (for instance, this can easily be shown with a tw-function of the form 
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\a — R\^^"' for large n). Besides, depending on the global surface density profile in the disc, it can be more 
appropriate to replace the homogeneous contribution (see Sect. 2) by an another analytical contribution 
(for instance E oc l/o instead of a constant as done here, provided kernels can be integrated). Finally, 
the method can easily be extended to treat non-axially symmetrical systems (Pierens & Hure 2005), and 
even tri-dimcnsional mass distributions. We plan to open soon to the scientific community an Internet site 
devoted to this question. 

The splitting method should be useful in many problems where great accuracy is nccdcid and/or fine 
physical effects must be investigated. For instance, we believe that significant progresses could be made in 
numerical simulations of planetary migrations which do show a great sensitivity to the field and potential 
values, even for low mass discs (e.g. Nelson & Benz 2003). 
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A. Splitting method for the potential 



According to notation, the general expression for the potential in the equatorial plane of a disc with 
inner edge ain and ttout is (e.g Durand, 1964): 



^ = -2G I 



;k^{a)K{k)da, 



(Al) 



For a constant surface density Sq, Eq.(Al) reads 



(•flout r~~ 

^homo.(^) = -2GSo / J-kK{k)da 

Jai^ V R 



(A2) 



For flin < R < flout, we separate this integral into two integrals, leftward and rightward to the fieldpoint 
where k — 1, namely 



homo. 



-2GSo 



%kK{k)da+ [ ,[^kK{k)da 

R Jr \ R 



(A3) 



Wc now set M = < 1 in the first integral and v = ^ < 1 in the second one. With the Gauss transformation 
(Gradshteyn & Ryzbik 1980) 

k(P^] ={l+u)K{u) u<l, 



Eq.(A2) becomes 



\l + uj 



/ uK{u)du~ I —^dv 

Jai„/R 



and so 



flout „l R \ 



E - E 

R V '^out / 



(if) 



(A4) 

(AS) 
(A6) 



B. Successive derivatives 



Prom Eq.(ll), we have 



da 
dw 



-2?i;aout > if R > a, 
if i? = a, 

2'u;aout > if R < a. 



(Bl) 



Let K be either or Kg, we have 



(ST, X K X \w\) = \w\ X {6Y, X k) ±ST, X K 
dw dw 

d da 
= \w\ X — (dS X k) X — ± dS X K 
da dw 

= 2w^aout X -;- (^5] x k) ± S'S x k 
da 



(B2) 



- 11 - 



for the first derivative, 



^(<5Ex.xH) = ^ 



2w^aout X -r (S'S X k) ±5T, X K 
da 



= Awttout-j^ {S^ X k) +4\wfal^i-^ {ST, x k)± 2|w|aout^ {S^ x k) 



for the second derivative is, and 



(B3) 



d d^ 
2 {2w ± \w\) aout^ (^S xk)+ 4|«;|3a2„,^ (5E x k) 



d d"^ 
6aout^ X k) + 8wal^^{\w\ ± 2w)^ (5S x k) 



d^ 
'da^ 



+ 8"'^«outi:^(^5]xK) 



(B4) 



for the third derivative. 



